% MotilityAnalysis_Displacement_Distance_Persistence.m
% Eden Ford
% 10 December 2023

% Code to visualize and analyze TrackMate data

% General definitions:
% Object - cell or cluster of cells that are touching/interacting with one
% another

close all
clear all
clc

%% User Inputs
FrameCutOff = 5; % anything 5 frames and less will not be included in the analysis

%% Data Import
% Distance traveled data
[num_0,txt_0,raw_0] = xlsread('MATLABCompiledData_0mM',1);
Distance0 = num_0;
[num_0,txt_0,raw_0] = xlsread('MATLABCompiledData_10mM',1);
Distance10 = num_0;

% Displacement data
[num_0,txt_0,raw_0] = xlsread('MATLABCompiledData_0mM',2);
Displacement0 = num_0;
[num_0,txt_0,raw_0] = xlsread('MATLABCompiledData_10mM',2);
Displacement10 = num_0;

%% Data Organization
% disp('TrackID  Time(hr)  Displacement(um)  DisplacementRate(um/hr)  Distance(um)  DistanceRate (um/hr)  Displacement:TotalDistance')
Ratio0 = [Distance0(:,1), Distance0(:,3), Displacement0(:,2), Displacement0(:,4), Distance0(:,2), Distance0(:,4), Displacement0(:,4)./Distance0(:,4)];
Ratio10 = [Distance10(:,1), Distance10(:,3), Displacement10(:,2), Displacement10(:,4), Distance10(:,2), Distance10(:,4), Displacement10(:,4)./Distance10(:,4)];

% Remove the stationary points
Ratio0_Clean = [];
Ratio10_Clean = [];
for i = 1:length(Ratio0)
    if isnan(Ratio0(i,7))
        continue
    elseif Ratio0(i,2) <= FrameCutOff/2
        continue
    else
        Ratio0_Clean = [Ratio0_Clean; Ratio0(i,:)];
    end
end
%Ratio0_Clean
for i = 1:length(Ratio10)
    if isnan(Ratio10(i,7))
        continue
    elseif Ratio10(i,2) <= FrameCutOff/2
        continue
    else
        Ratio10_Clean = [Ratio10_Clean; Ratio10(i,:)];
    end
end
% Ratio10_Clean
% DisplaceMax = [max(Ratio0_Clean(:,4)), max(Ratio10_Clean(:,4))]
% DistMax = [max(Ratio0_Clean(:,6)), max(Ratio10_Clean(:,6))]
% PersMax = [max(Ratio0_Clean(:,7)), max(Ratio10_Clean(:,7))]
%% Histograms

% Displacement Rate
figure
BinWidth = 5;
xMax = 40;
yMax = 1;
edges = linspace(0, 1, 21);
histogram_distribution_disp = zeros(2,20);

    % 0 mM CMP
    H1a = subplot(2,1,1);
    hold on
    histogram(Ratio0_Clean(:,4),'Normalization', 'probability', 'BinWidth', BinWidth);
    ylabel('Fraction of Objects')
    hold off
    title('0 mM CMP - Displacement Rate')
    xlim([0 xMax])
    ylim([0 yMax])
    histogram_distribution_disp(1,:) = histcounts(Ratio0_Clean(:,4),edges);

    % 10 mM CMP
    H1b = subplot(2,1,2);
    hold on
    histogram(Ratio10_Clean(:,4),'Normalization', 'probability', 'BinWidth', BinWidth);
    xlabel('Displacement Rate (um/hr)')
    ylabel('Fraction of Objects')
    hold off
    title('10 mM CMP - Displacement Rate')
    xlim([0 xMax])
    ylim([0 yMax])
    histogram_distribution_disp(2,:) = histcounts(Ratio10_Clean(:,4),edges);

% Total Distance Rate
figure
BinWidth = 5;
xMax = 60;
yMax = 0.6;
edges = linspace(0, 1, 21);
histogram_distribution_dist = zeros(2,20);

    % 0 mM CMP
    H2a = subplot(2,1,1);
    hold on
    histogram(Ratio0_Clean(:,6),'Normalization', 'probability', 'BinWidth', BinWidth);
    ylabel('Fraction of Objects')
    hold off
    title('0 mM CMP - Total Distance Rate')
    xlim([0 xMax])
    ylim([0 yMax])
    histogram_distribution_dist(1,:) = histcounts(Ratio0_Clean(:,6),edges);

    % 10 mM CMP
    H2b = subplot(2,1,2);
    hold on
    histogram(Ratio10_Clean(:,6),'Normalization', 'probability', 'BinWidth', BinWidth);
    xlabel('Total Distance Rate (um/hr)')
    ylabel('Fraction of Objects')
    hold off
    title('10 mM CMP - Total Distance Rate')
    xlim([0 xMax])
    ylim([0 yMax])
    histogram_distribution_dist(2,:) = histcounts(Ratio10_Clean(:,6),edges);

% Persistence
figure
BinWidth = 0.1;
xMax = 1;
yMax = 0.3;
edges = linspace(0, 1, 21);
histogram_distribution_per = zeros(2,20);

    % 0 mM CMP
    H3a = subplot(2,1,1);
    hold on
    histogram(Ratio0_Clean(:,7),'Normalization', 'probability', 'BinWidth', BinWidth);
    ylabel('Fraction of Objects')
    hold off
    title('0 mM CMP - Persistence')
    xlim([0 xMax])
    ylim([0 yMax])
    histogram_distribution_per(1,:) = histcounts(Ratio0_Clean(:,7),edges);

    % 10 mM CMP
    H3b = subplot(2,1,2);
    hold on
    histogram(Ratio10_Clean(:,7),'Normalization', 'probability', 'BinWidth', BinWidth);
    xlabel('Random Motility <--> Directional Persistence')
    ylabel('Fraction of Objects')
    hold off
    title('10 mM CMP - Persistence')
    xlim([0 xMax])
    ylim([0 yMax])
    histogram_distribution_per(2,:) = histcounts(Ratio10_Clean(:,7),edges);

%% Box Plot
% Displacement Rate
figure
DisplacementRate = [Ratio0_Clean(:,4); Ratio10_Clean(:,4)];
L1 = repmat({'0 mM mfCMP'},length(Ratio0_Clean),1);
L2 = repmat({'10 mM mfCMP'},length(Ratio10_Clean),1);
L = [L1; L2];
boxplot(DisplacementRate,L,'Notch','on')
xlabel('Hydrogel Composition')
ylabel('Displacement Rate (um/hr)')
title('Cell Displacement Rate')

% Total Distance Rate
figure
DistanceRate = [Ratio0_Clean(:,6); Ratio10_Clean(:,6)];
L1 = repmat({'0 mM mfCMP'},length(Ratio0_Clean),1);
L2 = repmat({'10 mM mfCMP'},length(Ratio10_Clean),1);
L = [L1; L2];
boxplot(DistanceRate,L,'Notch','on')
xlabel('Hydrogel Composition')
ylabel('Total Distance Rate (um/hr)')
title('Total Distance Rate')

% Persistence
figure
Persistence = [Ratio0_Clean(:,7); Ratio10_Clean(:,7)];
L1 = repmat({'0 mM mfCMP'},length(Ratio0_Clean),1);
L2 = repmat({'10 mM mfCMP'},length(Ratio10_Clean),1);
L = [L1; L2];
boxplot(Persistence,L,'Notch','on')
xlabel('Hydrogel Composition')
ylabel('Directional Persistence')
title('Directional Persistence')

%% Moments of Distribution
Type = {'0 mM mfCMP'; '10 mM mfCMP'};

% Displacement Rate
disp(' ');
disp('Cell Displacement Rate');
disp(' ');
ZeroCMP = Ratio0_Clean(:,4);
TenCMP = Ratio10_Clean(:,4);

AverageDisp = [mean(ZeroCMP); mean(TenCMP)]; % mean
StErrorDisp = [std(ZeroCMP)/sqrt(3); std(TenCMP)/sqrt(3)]; % standard error
MedianDisp = [median(ZeroCMP); median(TenCMP)]; % median
ModeDisp = [mode(ZeroCMP); mode(TenCMP)]; % mode
SkewnessDisp = [skewness(ZeroCMP); skewness(TenCMP)]; % skewness
KurtosisDisp = [kurtosis(ZeroCMP); kurtosis(TenCMP)]; % kurtosis

DisplacementTable = table(Type, AverageDisp, StErrorDisp, MedianDisp, ModeDisp, SkewnessDisp, KurtosisDisp);
DisplacementTable.Properties.VariableNames = {'GelComposition' 'Average' 'StError' 'Median' 'Mode' 'Skewness' 'Kurtosis'};
disp(DisplacementTable);

% Total Distance Rate
disp(' ');
disp('Total Distance Rate');
disp(' ');
ZeroCMP = Ratio0_Clean(:,6);
TenCMP = Ratio10_Clean(:,6);

AverageDist = [mean(ZeroCMP); mean(TenCMP)]; % mean
StErrorDist = [std(ZeroCMP)/sqrt(3); std(TenCMP)/sqrt(3)]; % standard error
MedianDist = [median(ZeroCMP); median(TenCMP)]; % median
ModeDist = [mode(ZeroCMP); mode(TenCMP)]; % mode
SkewnessDist = [skewness(ZeroCMP); skewness(TenCMP)]; % skewness
KurtosisDist = [kurtosis(ZeroCMP); kurtosis(TenCMP)]; % kurtosis

DistanceTable = table(Type, AverageDist, StErrorDist, MedianDist, ModeDist, SkewnessDist, KurtosisDist);
DistanceTable.Properties.VariableNames = {'GelComposition' 'Average' 'StError' 'Median' 'Mode' 'Skewness' 'Kurtosis'};
disp(DistanceTable);

% Persistence
disp(' ');
disp('Directional Persistence');
disp(' ');
ZeroCMP = Ratio0_Clean(:,7);
TenCMP = Ratio10_Clean(:,7);

AveragePers = [mean(ZeroCMP); mean(TenCMP)]; % mean
StErrorPers = [std(ZeroCMP)/sqrt(3); std(TenCMP)/sqrt(3)]; % standard error
MedianPers = [median(ZeroCMP); median(TenCMP)]; % median
ModePers = [mode(ZeroCMP); mode(TenCMP)]; % mode
SkewnessPers = [skewness(ZeroCMP); skewness(TenCMP)]; % skewness
KurtosisPers = [kurtosis(ZeroCMP); kurtosis(TenCMP)]; % kurtosis

PersistenceTable = table(Type, AveragePers, StErrorPers, MedianPers, ModePers, SkewnessPers, KurtosisPers);
PersistenceTable.Properties.VariableNames = {'GelComposition' 'Average' 'StError' 'Median' 'Mode' 'Skewness' 'Kurtosis'};
disp(PersistenceTable);
%% ANOVA
FillLength = length(Ratio10_Clean) - length(Ratio0_Clean);
Fill = NaN(FillLength,7);
Ratio0_NaN = [Ratio0_Clean; Fill];

% Displacement Rate
COL = 4;
DisplacementANOVA = [Ratio0_NaN(:,COL), Ratio10_Clean(:,COL)]; % matrix with data columns
[p_Disp tbl_Disp stats_Disp] = anova1(DisplacementANOVA)
DisplacementComparison = multcompare(stats_Disp) % multiple comparisons text

% Total Distance Rate
COL = 6;
DistanceANOVA = [Ratio0_NaN(:,COL), Ratio10_Clean(:,COL)]; % matrix with data columns
[p_Dist tbl_Dist stats_Dist] = anova1(DistanceANOVA)
DistanceComparison = multcompare(stats_Dist) % multiple comparisons text

% Persistence
COL = 7;
PersistenceANOVA = [Ratio0_NaN(:,COL), Ratio10_Clean(:,COL)]; % matrix with data columns
[p_Pers tbl_Pers stats_Pers] = anova1(PersistenceANOVA)
PersistenceComparison = multcompare(stats_Pers) % multiple comparisons text
